Table of contents

• Environment set up
• Dolomities area
• Nord-Ovest Italy Area
  • Region Level exploration
  • Municipality Level exploration

• CAI Alpine HUTs

  • Geocoding: from name to location
  • First exploration of the huts

• 4. Dolomities' Morphology

  • Explore the Catalogue Service Web
  • Lakes
  • Glaciers
  • Trentino Ortophoto

@Gabriele Ghisleni, Data Science.

---

  0. Environment set up

Install the required dependecies on google colab and mount the drive.

%%capture
!pip install geopandas==0.10.2
!pip install mapclassify==2.4.3
!pip install pygeos==0.10.2
!pip install rasterio==1.2.3
!pip install owslib==0.25.0
!pip install -U folium

import warnings, sys, os, folium, owslib, rasterio, geopy, shapely, json, requests, pprint
from folium import FeatureGroup, GeoJson, Marker
import matplotlib.pyplot as plt
from google.colab import drive
from tqdm.auto import tqdm
import geopandas as gpd
import pandas as pd
import numpy as np

drive.mount('/content/drive')
warnings.simplefilter("ignore")
p_print = pprint.PrettyPrinter(indent=4).pprint

# change this to your folder directory
folder_directory = "/content/drive/MyDrive/_____SHARED/Geospatial"
sys.path.append(folder_directory)

%%capture
#download the necessary utils files
if 'geo_utils' not in os.listdir():
  os.mkdir('geo_utils')

if "utils.py" not in os.listdir('geo_utils'):
  url_utils = "https://raw.githubusercontent.com/GabrieleGhisleni/Geo-Dolomities/master/utils/utils.py"
  url_style = "https://raw.githubusercontent.com/GabrieleGhisleni/Geo-Dolomities/master/utils/styles.py"
  !wget --no-cache --backups=1 -P'utils' {url_utils} 
  !wget --no-cache --backups=1 -P'utils' {url_style} 
  !touch utils/__init__.py
  assert 'utils.py' in os.listdir('geo_utils'), 'utils.py not present'
  assert 'styles.py' in os.listdir('geo_utils'), 'styles.py not present'

#in case ModuleError restart the kernel
import geo_utils.utils as utils 
import geo_utils.styles as styles
tiles = 'Stamen Terrain'

drive.mount('/content/drive')

folder_directory = "/content/drive/MyDrive/_____SHARED/Geospatial"
sys.path.append(folder_directory)
os.chdir(folder_directory)

import utils, styles
tiles = 'Stamen Terrain'

Mounted at /content/drive

  1. Dolomities area

Let's explore the area of the dolomities, the data are from the offical website of Dolomiti Unesco, that are available in .kml format.

We load the data and perform some preprocessing before so to have a better insight.

# read kml in GeoPandas https://docs.astraea.earth/hc/en-us/articles/360043923831-Read-a-KML-File-into-a-GeoPandas-DataFrame

gpd.io.file.fiona.drvsupport.supported_drivers['KML'] = 'rw'
dolomiti_df = gpd.read_file('data/geo_dolomiti/dolomiti_area.kml', driver='KML')

print(f"Reference system information:\n{ '-'*50}\n{dolomiti_df.crs.__repr__()}{'-'*50}")

Reference system information:
--------------------------------------------------
<Geographic 2D CRS: EPSG:4326>
Name: WGS 84
Axis Info [ellipsoidal]:
- Lat[north]: Geodetic latitude (degree)
- Lon[east]: Geodetic longitude (degree)
Area of Use:
- name: World.
- bounds: (-180.0, -90.0, 180.0, 90.0)
Datum: World Geodetic System 1984 ensemble
- Ellipsoid: WGS 84
- Prime Meridian: Greenwich
--------------------------------------------------

# some preprocessing
dolomiti_df["name"] = dolomiti_df["Name"].apply(lambda x: x.replace(',', ' &').replace('\xa0', ' '))
dolomiti_df.loc[2, 'name'] = "Sistema 3 - Pale di San Martino"
dolomiti_df['url_info'] = dolomiti_df.Description.apply(lambda x: x.split('>')[1])

# representative points and area expressed in m^2 (project to 32632 which is in meters).
dolomiti_df['representative_point'] = dolomiti_df.to_crs(epsg=32632).geometry.representative_point().to_crs(epsg=4326)
dolomiti_df['representative_point'] = dolomiti_df['representative_point'].apply(str)
dolomiti_df['area'] = dolomiti_df.to_crs(epsg=32632).geometry.area.apply(lambda x: int(x)) #metersquare since we are in 32632

dolomiti_df = dolomiti_df.loc[:, ['name', 'url_info', 'area', 'representative_point', 'geometry']]
dolomiti_df.head(3)

	name	url_info	area	representative_point	geometry
0	Sistema 1 - Pelmo & Croda da Lago	http://www.dolomitiunesco.info/?gruppo-dolomit...	43153031	POINT (12.08268734287994 46.45014333277591)	POLYGON Z ((12.12142 46.40488 0.00000, 12.1214...
1	Sistema 2 - Marmolada	http://www.dolomitiunesco.info/?gruppo-dolomit...	22057860	POINT (11.859039731003216 46.428646792973986)	POLYGON Z ((11.82319 46.45296 0.00000, 11.8221...
2	Sistema 3 - Pale di San Martino	http://www.dolomitiunesco.info/?gruppo-dolomit...	316567741	POINT (12.165681166978791 46.23429736260587)	POLYGON Z ((11.98683 46.31720 0.00000, 11.9880...

# Some information regarding the groups and their extensions
for idx, row in dolomiti_df.sort_values('area', ascending=False).reset_index().iterrows():
  print(f"""\
{idx+1}° largest group: '{row["name"]}' Area: {round(row.area * 1e-6)} km2.
Website for more info: {row.url_info}.
  """)

1° largest group: 'Sistema 5 - Dolomiti settentrionali' Area: 536 km2.
Website for more info: http://www.dolomitiunesco.info/?gruppo-dolomitico=dolomiti-settentrionali.
  
2° largest group: 'Sistema 3 - Pale di San Martino' Area: 317 km2.
Website for more info: http://www.dolomitiunesco.info/?gruppo-dolomitico=pale-di-san-martino.
  
3° largest group: 'Sistema 4 - Dolomiti friulane e d'Oltre Piave' Area: 215 km2.
Website for more info: http://www.dolomitiunesco.info/?gruppo-dolomitico=dolomiti-friulane.
  
4° largest group: 'Sistema 9 - Dolomiti di Brenta' Area: 111 km2.
Website for more info: http://www.dolomitiunesco.info/?gruppo-dolomitico=dolomiti-di-brenta.
  
5° largest group: 'Sistema 7 - Sciliar-Catinaccio & Latemar' Area: 93 km2.
Website for more info: http://www.dolomitiunesco.info/?gruppo-dolomitico=sciliar-catinaccio-latemar.
  
6° largest group: 'Sistema 6 - Puez-Odle' Area: 79 km2.
Website for more info: http://www.dolomitiunesco.info/?gruppo-dolomitico=puez-odle.
  
7° largest group: 'Sistema 1 - Pelmo & Croda da Lago' Area: 43 km2.
Website for more info: http://www.dolomitiunesco.info/?gruppo-dolomitico=pelmo-croda-da-lago.
  
8° largest group: 'Sistema 2 - Marmolada' Area: 22 km2.
Website for more info: http://www.dolomitiunesco.info/?gruppo-dolomitico=marmolada.
  
9° largest group: 'Sistema 8 - Bletterbach' Area: 3 km2.
Website for more info: http://www.dolomitiunesco.info/?gruppo-dolomitico=bletterbach.
  

# first representation and plot of the dolomities
fig, axes = plt.subplots(1,1, figsize=(15,15))

dolomiti_df.to_crs(epsg=4326).plot(
    ax=axes,
    column='name',
    cmap = "tab10",
    legend=True,  
    categorical=True,
    legend_kwds= {
        'bbox_to_anchor':(.3, 1.05),
        'fontsize':15,
        'frameon':False
        }
)

plt.title('Dolomities Area', fontdict=dict(size=25))
plt.tight_layout()
plt.axis('off')
plt.show()

  2. Nord-Ovest Italy Area

Extract the regions of the italy that are closeby the dolomities.

  2.1 Region level exploration

italy_regions = gpd.read_file("data/istat_administrative.gpkg", layer="regions")

# Set correct crs
italy_regions = italy_regions.to_crs("EPSG:4326")

# Selecting only interested regions
assert italy_regions.crs == dolomiti_df.crs

selected_regions = []

for _, row in dolomiti_df.iterrows():
  regions = italy_regions.geometry.overlaps(row.geometry)
  selected_regions.extend(
      italy_regions.loc[regions].DEN_REG.values.tolist()
  )
  
selected_regions = list(set(selected_regions))
dolomiti_regions = italy_regions.loc[italy_regions.DEN_REG.isin(selected_regions)]

Using Folium library for creating nice interactive maps so to explore the dolomities with respect to the italian regions.

# https://gis.stackexchange.com/questions/392531/modify-geojson-tooltip-format-when-using-folium-to-produce-a-map
# https://stackoverflow.com/questions/58305337/plotting-polygons-in-python-using-geopandas-and-folium
# https://stackoverflow.com/questions/63960271/folium-geojson-custom-color-map
# https://leafletjs.com/reference.html#path-option

dolomiti_map = utils.get_new_map(title='Dolomities in Italy', tiles=tiles)

feature_layers = dict(
  dolomiti_area = FeatureGroup(name='Area Dolomiti'),
  regions_layer = FeatureGroup(name='Italian Regions')   
)

dolomiti_region = GeoJson(
    data = dolomiti_regions,
    style_function = styles.style_region,
    tooltip = folium.GeoJsonTooltip(fields=['DEN_REG'])
  ).add_to(feature_layers["regions_layer"])

dolomiti_geo = GeoJson(
    data = dolomiti_df,
    style_function = styles.style_dolomiti,
    tooltip = folium.GeoJsonTooltip(
            fields=[
                'name',
                'url_info',
                'area',
                'representative_point'
            ])
  )

dolomiti_geo.add_to(feature_layers["dolomiti_area"])
utils.add_map_infos(feature_layers, dolomiti_map, reverse_order=True)

Make this Notebook Trusted to load map: File -> Trust Notebook

# explore with static maps
fig, axes = plt.subplots(1,1, figsize=(15,15))

dolomiti_regions.to_crs(epsg=4326).plot(
    ax=axes,
    cmap = "Pastel1",
    column='DEN_REG',
    categorical=True,
    legend=True,  
  )

dolomiti_df.to_crs(epsg=4326).plot(
    ax=axes,
    column='name',
    cmap = "tab10",
    legend=True,  
    categorical=True,
    legend_kwds= {
        'bbox_to_anchor':(.3, 1.05),
        'fontsize':10,
        'frameon':False
        }
  )

plt.title('Dolomities Area', fontdict=dict(size=25))
plt.tight_layout()
plt.axis('off')
plt.show()

  2.2 Muncipality level exploration

italy_municipalities = gpd.read_file("data/istat_administrative.gpkg", layer="municipalities")

# filter by the interested regions
italy_municipalities = italy_municipalities.loc[italy_municipalities.COD_REG.isin(dolomiti_regions.COD_REG.unique())]

# Set correct crs
italy_municipalities = italy_municipalities.to_crs("EPSG:4326")

# find the municipalities that are nearby the dolomities area
def geometry_in_areas(geom, regions_df, attribute_name='name'):
  res = []
  for idx, region in regions_df.iterrows():
    geometry = region.geometry
    if geom.overlaps(geometry): 
      res.append(region[attribute_name])

  return ','.join(res)
  
buffered_dolomiti = dolomiti_df.copy()
buffered_dolomiti.geometry = dolomiti_df.to_crs(epsg=32632).geometry.buffer(500).to_crs(epsg=4326).values

selected_municipalities = italy_municipalities.geometry.apply(lambda x: geometry_in_areas(x, buffered_dolomiti)) # dolomiti_df dolomiti_regions
italy_municipalities['dolomiti'] = selected_municipalities.values
dolomiti_municipalities = italy_municipalities.loc[selected_municipalities != '']
dolomiti_municipalities = dolomiti_municipalities.loc[:, ['COD_REG', 'COMUNE', 'geometry', 'dolomiti']]

dolomiti_municipalities.head(2)

	COD_REG	COMUNE	geometry	dolomiti
2760	4	Badia	MULTIPOLYGON (((11.89993 46.63847, 11.89978 46...	Sistema 5 - Dolomiti settentrionali,Sistema 6 ...
2763	4	Braies	MULTIPOLYGON (((12.13557 46.73782, 12.13594 46...	Sistema 5 - Dolomiti settentrionali

# add regions of the municipality
dolomiti_municipalities = dolomiti_municipalities.merge(
    dolomiti_regions.loc[:, ['COD_REG', 'DEN_REG']],
    on='COD_REG'
  )

# municipality for each region that touches the dolomities
for _, row in dolomiti_municipalities.groupby('DEN_REG').size().sort_values(0, ascending=False).reset_index().iterrows():
  print(f"'{row.DEN_REG}' has {row[0]} municipalities nearby the dolomities.")

'Trentino-Alto Adige' has 44 municipalities nearby the dolomities.
'Veneto' has 40 municipalities nearby the dolomities.
'Friuli Venezia Giulia' has 6 municipalities nearby the dolomities.

# create the map and add geojson infos
dolomiti_map_municipalities = utils.get_new_map(title='Dolomities over municipalities', tiles='OpenStreetMap')

feature_layers = dict(
  dolomiti_area = FeatureGroup(name='Area Dolomiti'),
  municipalities_layer = FeatureGroup(name='Italian Municipalities')
)

dolomiti_geo.add_to(feature_layers["dolomiti_area"])

GeoJson(
    data = dolomiti_municipalities,
    style_function = styles.style_municipalities,
    tooltip = folium.GeoJsonTooltip(
            fields=[
                'COMUNE',
                'DEN_REG',
                'dolomiti'
            ])
).add_to(feature_layers["municipalities_layer"])

utils.add_map_infos(feature_layers, dolomiti_map_municipalities,  reverse_order=True)

Make this Notebook Trusted to load map: File -> Trust Notebook

  3. CAI Alpine HUT

Manually extract the CAI from the CAI official website alpine hut and organize them into a useful way.

# https://www.cai.it/andare-in-montagna/rifugi-e-bivacchi/rifugi-del-club-alpino-italiano/

rifugi_trentino = pd.read_excel(
    io='data/alpine_huts/rifugi_trentino.xlsx', 
    sheet_name=0, 
    header=1
  )

rifugi_veneto = pd.read_excel(
    io='data/alpine_huts/rifugi_veneto.xlsx', 
    sheet_name=0, 
    header=1
  )

rifugi = pd.concat([rifugi_trentino, rifugi_veneto])
rifugi.head(2)

	Nome	Tipo	Sezione	Quota (m slm)	Posti totali	Categoria
0	Rifugio Vioz	Rifugio custodito	S.A.T.	3535	66	E
1	Rifugio Boè	Rifugio custodito	S.A.T.	2873	69	D

  3.1 Geocoding: from name to location

user_agent = "Mozilla/5.0 (Linux; Android 10) AppleWebKit/537.36 (KHTML, like Gecko) Chrome/86.0.4240.75 Mobile Safari/537.36"
geolocator = geopy.geocoders.Nominatim(user_agent=user_agent, timeout=3)

def geocode_sat(rifugi):
  detail = []
  not_founded = []

  for idx, row in tqdm(rifugi.iterrows(), total=len(rifugi)):
    name = row['Nome']
    # manually fix the names that differ
    if name in utils.names_adj:
      name = utils.names_adj[name]

    # geolocate the huts
    geocodes = geolocator.geocode(
            query = name, 
            addressdetails = True,
            extratags = True
          )
    
    if geocodes:
      detail.append(geocodes.raw)
    else:
      not_founded.append(row['Nome'])

  rifugi = utils.expand_raw_geocode(pd.DataFrame(detail))
  rifugi = gpd.GeoDataFrame(rifugi, crs="EPSG:4326")
  rifugi = rifugi.loc[~rifugi.county.isin(['Roma Capitale', "Reggio nell'Emilia", 'Bergamo'])]
  rifugi.drop(columns=['boundingbox']).to_file('data/rifugi.shp') 
  return rifugi
  

cai_rifugi = (
    gpd.read_file("data/alpine_huts/rifugi.shp") 
    if os.path.exists('data/alpine_huts/rifugi.shp')
    else geocode_sat(rifugi)
  )

cai_rifugi.head(2)

	municipali	county	state	road_name	extended_n	opening_ho	contact_mo	winter_roo	capacity	website	email	display_na	type_descr	geometry
0	Comunità della Valle di Sole	Provincia di Trento	Trentino-Alto Adige/Südtirol	Vioz	Rifugio Mantova al Vioz	Jun 20-Sep 20	+39 339 2798826	yes	30	http://www.rifugiovioz.it/	info@rifugiovioz.it	Rifugio Mantova al Vioz, Vioz, Peio, Comunità ...	alpine_hut	POINT (10.63578 46.39920)
1	Pustertal - Val Pusteria	Bolzano - Bozen	Trentino-Alto Adige/Südtirol	Cresta Strenta - Lichtenfelser Weg	Rifugio Boé	Jun 20-Sep 20	+39 349 5730110	yes	69	http://www.rifugioboe.it/	rifugioboe@gmail.com	Rifugio Boé, Cresta Strenta - Lichtenfelser We...	alpine_hut	POINT (11.82328 46.51461)

  3.2 First exploration of the huts

sat_map = utils.get_new_map(title='Hut mantained by SAT', tiles=tiles)

unique_provice_dolomities = cai_rifugi.county.unique().tolist()
colors = ['gray', 'lightgreen', 'green', 'orange', 'darkred','black']

feature_layers = {
    group:FeatureGroup(name=f"<span style='color:{utils.c_colors[idx]}'>Hut in {group}</span>") 
    for idx, group in enumerate(unique_provice_dolomities)
  }

feature_layers["dolomiti_area"] = FeatureGroup(name='Area Dolomiti')
dolomiti_geo.add_to(feature_layers["dolomiti_area"])

for idx, row in cai_rifugi.iterrows():
  marker = Marker(
      location = [i[0] for i in reversed(row.geometry.xy)],
      popup = utils.get_info_from_row(row),
      tooltip = utils.get_info_from_row(row),
      icon = folium.Icon(color= utils.c_colors[unique_provice_dolomities.index(row.county)])
  )

  marker.add_to(feature_layers[row.county])


utils.add_map_infos(feature_layers, sat_map)

Make this Notebook Trusted to load map: File -> Trust Notebook

# explore with static maps
fig, axes = plt.subplots(1,1, figsize=(15,15))

dolomiti_regions.to_crs(epsg=4326).plot(ax=axes,  linewidth=3, facecolor='none', edgecolor='gray')
dolomiti_df.to_crs(epsg=4326).plot(
    ax=axes,
    linewidth=2,
    edgecolor='black',
    column='name',
    cmap = "tab10",
    legend=True,  
    categorical=True,
    legend_kwds= {
        'bbox_to_anchor':(.3, 1.05),
        'fontsize':10,
        'frameon':False
        }
)

cai_rifugi.to_crs(epsg=4326).plot(ax=axes, color = "orange")

plt.title('Huts over the Dolomities area', fontdict=dict(size=25))
plt.tight_layout()
plt.axis('off')
plt.show()

  4. Dolomities' Morphology

  4.1 Explore the Catalogue Service Web

from owslib.csw import CatalogueServiceWeb

# fetch the catalog
csw = CatalogueServiceWeb("http://www.pcn.minambiente.it/geoportal/csw")
query = owslib.fes.PropertyIsLike('csw:AnyText', 'Laghi')
csw.getrecords2(constraints=[query], maxrecords=1_000)

#explore the available data
for i, rec in enumerate(csw.records): 
  print(f"{rec} - {csw.records[rec].title}")
  if i > 15: break

r_campan:16050024:20090806:141714 - Specchi acqua campania
r_piemon:4ba5287d-0c81-4729-84b9-906ad5ccb173 - Idrografia 1:100.000 - Laghi  - Storico
r_vda:03151-META:20170302:081500 - Arpa - Catasto laghi della Valle d'Aosta
r_piemon:60048c12-904a-4aa5-a062-1a9dd94ae971 - Ppr - Laghi (tavv. P4-P5)
p_TN:0403efca-cadb-4463-a656-4d0fba207953 - Aree di protezione dei laghi - Variante PUP 2000
r_piemon:4fb63c5b-813b-4c6d-9fed-6ce92df6ca64 - Ppr - Laghi (tav. P2)
r_basili:54eecefc:162a68d9be7:-2723 - Beni paesaggistici art. 142 let. b - Laghi ed invasi artificiali (Buffer)
ispra_rm:20101021:113000 - Reticolo Idrografico Nazionale - Dataset
r_lombar:36BB78E0-8946-4E7D-ADAC-7D149D18B570 - PTUA 2016 - Tav. 1 - Corpi idrici superficiali e bacini drenanti (fiumi e laghi)
R_SARDEG:XBDOJ - PPR06 - Laghi naturali, invasi artificiali, stagni e lagune
r_campan:16050003:20090722:144054 - Aree di tutela coste e laghi
r_piemon:b65fc959-dfa5-4edd-991d-b6b9f638ab5b - Ppr - Lettera b) fasce lacuali (tav. P2)
p_TN:9a5dc19c-543f-44d3-a36c-b4ea73d0a48a - Z101_P_PUP
abnfa_tn:00002:20120504:130718 - Laghi Alpi Orientali
p_TN:822c7446-8e84-45b6-9d50-98d5dca19416 - Z310_P_PUP
r_sicili:aa3781a2-3c7d-4a32-ac9d-e8aae4195f1f - Carta dell' uso dei suoli  1:100.000
r_piemon:8d267850-e543-4294-bb23-483aa5335183 - Corpi idrici WFD 2000/60/CE Triennio 2009-2011 - Geo servizio WFS

  4.2 Lakes

Explore the lakes in the dolomities.

# explore the resource found
unique_id = "m_amte:299FN3:5237cf1e-4d49-486c-a688-152bc4473508"
dicts = csw.records[unique_id].__dict__

print('Abstract:\n'); print(dicts['abstract']); print('-'*100)
print('References:\n'); p_print(dicts['references']); print('-'*100)
print('Information:\n'); p_print({f"{k}: {v}" for k,v in dicts.items() if v and k != 'xml' and k != 'abstract' and type(v) != list}); print('-'*100)

Abstract:

Elementi poligonali appartenenti al DBPrior10K  rappresentanti gli specchi dacqua interni come laghi, lagune, stagni, bacini artificiali, ecc.
----------------------------------------------------------------------------------------------------
References:

----------------------------------------------------------------------------------------------------
Information:

----------------------------------------------------------------------------------------------------

from owslib.wfs import WebFeatureService

url = "http://wms.pcn.minambiente.it/ogc?map=/ms_ogc/WMS_v1.3/Vettoriali/Specchi_Acqua.map"
wfs = WebFeatureService(url=url, version="1.1.0")

{f"{k}: {v}" for k,v in wfs.__dict__.items() if v and k != 'xml' and k != 'abstract' and type(v) != list}

{'_capabilities: <Element {http://www.opengis.net/wfs}WFS_Capabilities at 0x7f061f7c6460>',
 'auth: <Authentication shared=False username=None password=None cert=None verify=True auth_delegate=None>',
 "contents: {'ID.ACQUEFISICHE.SPECCHI.ACQUA': <owslib.feature.wfs110.ContentMetadata object at 0x7f061e43e9d0>}",
 'filters: <owslib.fes.FilterCapabilities object at 0x7f061d908b90>',
 'identification: Service: OGC WFS, title=',
 'owscommon: <owslib.ows.OwsCommon object at 0x7f061d90e250>',
 'provider: <owslib.ows.ServiceProvider object at 0x7f061d9846d0>',
 'timeout: 30',
 'url: http://wms.pcn.minambiente.it/ogc?map=%2Fms_ogc%2FWMS_v1.3%2FVettoriali%2FSpecchi_Acqua.map',
 'version: 1.1.0'}

# load the layer as geodataframe

lakes_df = gpd.read_file(
    wfs.getfeature(
    typename="ID.ACQUEFISICHE.SPECCHI.ACQUA",
    bbox = dolomiti_df.total_bounds.tolist(),
  )
)


# invert the xy coordinates
def reverse_geom(geom):
    def _reverse(x, y):
        return y[::-1], x[::-1]

    return shapely.ops.transform(_reverse, geom)


lakes_df.geometry = lakes_df.geometry.apply(reverse_geom)

buffered_dolomiti = dolomiti_df.copy()
buffered_dolomiti.geometry = dolomiti_df.to_crs(epsg=32632).geometry.buffer(5_000).to_crs(epsg=4326).values

# find the lakes in the dolomities area
idxs = []
for idx, row in lakes_df.to_crs(epsg=4326).iterrows():
  if buffered_dolomiti.geometry.contains(row.geometry).sum() > 0:
    idxs.append(idx)

lakes_df = lakes_df.iloc[idxs]
lakes_df = lakes_df.loc[lakes_df.nome != ""].reset_index(drop=True)
lakes_df["area"] = lakes_df.to_crs(epsg=32632).geometry.area
lakes_df["area"] = lakes_df.to_crs(epsg=32632).geometry.area
lakes_df["area"] = (lakes_df["area"] * 1e-6).apply(lambda x: round(x, 3))


#biggest lakes
for _, row in lakes_df.sort_values('area', ascending=False).head().iterrows():
  print(f"Name: {row['nome']} area: {row['area']}")

Name: LAGO DI MOLVENO area: 3.218
Name: Lago del Mis area: 1.098
Name: LAGO ARTIFICIALE DI FEDAIA area: 0.545
Name: LAGO DI TOVEL area: 0.371
Name: LAGO DELLO SCHENER area: 0.356

  4.3 Glaciers

Exploration of all the glaciers that are in the dolomities.

data from http://www.glaciologia.it/en/, thanks to:

Salvatore M.C., Zanoner T., Baroni C., Carton A., Banchieri F.A., Viani C., Giardino M., Perotti L. (2015) - The state of Italian glaciers: A snapshot of the 2006-2007 hydrological period. Geografia Fisica e Dinamica Quaternaria, 38 (2), pp. 175-198.

# load the data if present otherwhise download 
if 'ghiacciai.json' not in os.listdir('data'):
  raw_glacier = requests.get('https://repo2.igg.cnr.it/ghiacciaiCGI/data/datagralp.json').json() 

else:
  with open('data/ghiacciai.json') as f:
    raw_glacier = json.load(f)

# set the correct crs
crs = raw_glacier['crs']['properties']['name'].split(':')[-1]
glacier_df =  gpd.GeoDataFrame.from_features(raw_glacier, crs=crs)

# keeps only the glaciers into the dolomities
idxs = []
for idx, row in glacier_df.to_crs(epsg=4326).iterrows():
  if dolomiti_df.geometry.contains(row.geometry).sum() > 0:
    idxs.append(idx)

glacier_df = glacier_df.iloc[idxs].reset_index(drop=True)
glacier_df["area"] = glacier_df.to_crs(epsg=32632).geometry.area
glacier_df["area"] = (glacier_df["area"] * 1e-6).apply(lambda x: round(x, 3))


#biggest glacier
for _, row in glacier_df.sort_values('area', ascending=False).head().iterrows():
  print(f"Name: {row['GLACIER']} area: {row['area']}, gruppo: {row['GROUP_NAME']}")

Name: Ghiacciaio Principale della Marmolada area: 1.48, gruppo: Marmolada
Name: Ghiacciaio Superiore (Orientale) dell'Antelao area: 0.255, gruppo: Antelao
Name: Ghiacciaio (Vedretta) del Lagol o de l'Agola area: 0.215, gruppo: Brenta
Name: Ghiacciaio di Travignolo area: 0.194, gruppo: Pale di San Martino
Name: Ghiacciaio Occidentale del Sorapis area: 0.179, gruppo: Sorapiss

areas_glacier = glacier_df.groupby('GROUP_NAME').sum().reset_index().loc[:, ['GROUP_NAME', 'area']].sort_values('area', ascending=False)
for _, row in areas_glacier.iterrows():
  print(f'Glacier group: {row.GROUP_NAME} has {round(row.area, 3)} km^2 of glaciers ')

Glacier group: Marmolada has 2.025 km^2 of glaciers 
Glacier group: Brenta has 1.124 km^2 of glaciers 
Glacier group: Pale di San Martino has 0.602 km^2 of glaciers 
Glacier group: Sorapiss has 0.453 km^2 of glaciers 
Glacier group: Antelao has 0.404 km^2 of glaciers 
Glacier group: Marmarole has 0.384 km^2 of glaciers 
Glacier group: Cristallo has 0.335 km^2 of glaciers 
Glacier group: Popera has 0.263 km^2 of glaciers 
Glacier group: Pelmo has 0.176 km^2 of glaciers 
Glacier group: Civetta has 0.109 km^2 of glaciers 
Glacier group: Tofane has 0.062 km^2 of glaciers 
Glacier group: Fanis has 0.037 km^2 of glaciers 
Glacier group: Croda Rossa d'Ampezzo has 0.001 km^2 of glaciers 
Glacier group: Moiazza has 0.001 km^2 of glaciers 

# explore the lakes and the glaciers in the dolomities
glacier_lake_map = utils.get_new_map(title='Glaciers and Lakes', tiles='Stamen Terrain')

feature_layers = dict(
  dolomiti_area = FeatureGroup(name='Area Dolomiti'),
  glacier = FeatureGroup(name='Glacier'),
  lakes = FeatureGroup(name='Lakes'),
)

dolomiti_geo.add_to(feature_layers["dolomiti_area"])

GeoJson(
    data = glacier_df,
    style_function = lambda x : {'color': 'purple', 'fill': True},
    tooltip = folium.GeoJsonTooltip(
              fields=[
                  'GLACIER',
                  'GROUP_NAME',
                  'MAX_ELEV',
                  'MIN_ELEV',
                  'DATE',
                  'AREA_KM'
              ])
  ).add_to(feature_layers["glacier"])


GeoJson(
    data = lakes_df,
    style_function = lambda x : {'color': 'blue', 'fill': True},
    tooltip = folium.GeoJsonTooltip(
              fields=[
                  'nome',
                  'area'
              ])
  ).add_to(feature_layers["lakes"])

utils.add_map_infos(feature_layers, glacier_lake_map)

Make this Notebook Trusted to load map: File -> Trust Notebook

# explore the lakes and the glaciers in the dolomities
fig, axes = plt.subplots(1,1, figsize=(15,15))

dolomiti_df.to_crs(epsg=4326).plot(
    ax=axes,
    linewidth=2,
    edgecolor='black',
    column='name',
    cmap = "tab10",
    alpha=.15,
    legend=True,  
    categorical=True,
    legend_kwds= {
        'bbox_to_anchor':(.3, 1.05),
        'fontsize':10,
        'frameon':False
        }
)

glacier_df.to_crs(epsg=4326).plot(ax=axes, color = "purple")
lakes_df.to_crs(epsg=4326).plot(ax=axes, color = "blue")

plt.title('Glaciers and Lakes in the Dolomities area', fontdict=dict(size=25))
plt.tight_layout()
plt.axis('off')
plt.show()

  4.4 Trentino Ortophoto

Quickly exploration of the ortophoto of trentino only.

La Provincia Autonoma di Trento ha la piena proprietà dell'ortofoto. Attribuzione 3.0 Italia (CC BY 3.0 IT).

# fetch the catalog
csw = CatalogueServiceWeb("http://geodati.gov.it/RNDT/csw")
trento_query = owslib.fes.PropertyIsLike('csw:AnyText', 'Trento')
csw.getrecords2(constraints=[trento_query], maxrecords=1_000)

#explore the available data
for rec in csw.records: 
  if 'orto' in csw.records[rec].title.lower():
    print(f"{rec} - {csw.records[rec].title}")

agea:00129:20090724:090446 - Ortofotocarta Trento 2003
agea:00377:20090911:093144 - Ortofotocarta Trento 2008
agea:00128:20090724:085449 - Ortofotocarta Trento 1997
c_l378:ortofoto2009 - Ortofoto 2009
p_TN:fbbc1e07-0b8e-46c9-b961-a02d8bebb217 - Ortofoto PAT 1973 in scala di grigi
p_TN:013ef530-ee77-49d2-8f95-035b27ab1f0a - Ortofoto Val di Sella 2007 RGB
p_TN:44cd8a88-ba6e-4299-8075-df9996401006 - Ortofoto Padergnone 2007 RGB
p_TN:e506ab2b-5eb6-452a-8ebb-102ffc4454b1 - Ortofoto Paneveggio 2007 RGB
p_TN:054e1af3-2dd8-4496-8b83-261cccc2674c - Ortofoto Ravina 2007 RGB
p_TN:f2e88f1b-05d9-4942-93ee-857a0a9e1f0b - Ortofoto PAT 2015 RGB

from owslib.wms import WebMapService
from rasterio.plot import show

orto_photo_path = "p_TN:f2e88f1b-05d9-4942-93ee-857a0a9e1f0b"
wms_trento = csw.records[orto_photo_path].__dict__['references'][3]['url']
wms = WebMapService(wms_trento)

trentino_map = wms.getmap(
    layers=['ecw-rgb-2015'],
    srs='EPSG:25832',
    format='image/jpeg',
    bbox=(609500.0, 5056790.0, 731610.0, 5158990.0),
    size=(3500,3500)
  )

fig, axes = plt.subplots(1,1, figsize=(15,15))
fig = show(
    rasterio.MemoryFile(trentino_map).open(),
    ax=axes, 
    title='Trentino Natural View'
)

plt.xticks([]); plt.yticks([])
plt.show()